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Abstract 


The need for accurate material models to simulate the deformation, damage and failure of polymer 
matrix composites under impact conditions is becoming critical as these materials are gaining increased 
use in the aerospace and automotive communities. The aerospace community has identified several key 
capabilities which are currently lacking in the available material models in commercial transient dynamic 
finite element codes. To attempt to improve the predictive capability of composite impact simulations, a 
next generation material model is being developed for incorporation within the commercial transient 
dynamic finite element code LS-DYNA. The material model, which incorporates plasticity, damage and 
failure, utilizes experimentally based tabulated input to define the evolution of plasticity and damage and 
the initiation of failure as opposed to specifying discrete input parameters such as modulus and strength. 
The plasticity portion of the orthotropic, three-dimensional, macroscopic composite constitutive model is 
based on an extension of the Tsai-Wu composite failure model into a generalized yield function with a 
nonassociative flow rule. For the damage model, a strain equivalent formulation is used to allow for the 
uncoupling of the deformation and damage analyses. In the damage model, a semicoupled approach is 
employed where the overall damage in a particular coordinate direction is assumed to be a multiplicative 
combination of the damage in that direction resulting from the applied loads in various coordinate 
directions. For the failure model, a tabulated approach is utilized in which a stress or strain based 
invariant is defined as a function of the location of the current stress state in stress space to define the 
initiation of failure. Failure surfaces can be defined with any arbitrary shape, unlike traditional failure 
models where the mathematical functions used to define the failure surface impose a specific shape on the 
failure surface. In the current paper, the complete development of the failure model is described and the 
generation of a tabulated failure surface for a representative composite material is discussed. 


Introduction 


As composite materials are gaining increased use in aircraft components where impact resistance is 
critical (such as the turbine engine fan case), the need for accurate material models to simulate the 
deformation, damage and failure response of polymer matrix composites under impact conditions is 
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gaining in importance. While there are several material models currently available within commercial 
transient dynamic codes such as LS-DYNA (Ref. 1) to analyze the impact response of composites, areas 
have been identified where the predictive capabilities of these models can be improved. Most importantly, 
the existing models often require extensive correlation based on structural level impact tests, which 
significantly limits the ability to use these models as predictive tools. Furthermore, most of the existing 
models apply either a plasticity based approach (such as that used by Sun and Chen (Ref. 2)) or a 
continuum damage mechanics approach (such as that used by Matzenmiller et al. (Ref. 3)) to simulate the 
nonlinearity that takes place in the composite response. As documented in detail by Goldberg, et al. 
(Refs. 4 and 5), either of these approaches can capture certain aspects of the actual composite behavior. 
However, optimally, combining a plasticity based deformation model (to capture the rate dependence and 
significant nonlinearity, particularly in shear, observed in the composite response) with a damage model 
(to account for the changes in the unloading modulus observed as the material is unloaded from various 
stress levels, as well as to account for strain softening observed after the peak stress is reached) can 
provide advantages over using one approach or the other. In addition, the input to current material models 
generally consists of point-wise properties (such as the modulus, failure stress or failure strain in a 
particular coordinate direction) that leads to curve fit approximations to the material stress-strain curves. 
This type of approach either leads to models with only a few parameters, which provide a crude 
approximation at best to the actual material response, or to models with many parameters which require a 
large number of complex tests to characterize. An improved approach is to use tabulated data, obtained 
from a well-defined set of straightforward experiments. Using tabulated data allows the actual material 
response data to be entered in a discretized form, which permits a more accurate representation of the 
actual material response. 

To begin to address these needs, a new composite material model is being developed and 
implemented for use within LS-DYNA. The material model is meant to be a fully generalized model 
suitable for use with any composite architecture (laminated or textile). The deformation model is based on 
extending the commonly used Tsai-Wu composite failure model (Ref. 6) to a strain hardening plasticity 
model with a nonassociative flow rule. For the damage model, a strain equivalent formulation is used in 
which the deformation and damage calculations can be uncoupled. A significant feature in the developed 
damage model is that a semicoupled approach has been utilized in which a load in a particular coordinate 
direction results in damage (and thus stiffness reduction) in multiple coordinate directions. While 
different from the approach used in many existing damage mechanics models (Ref. 3), in which a load in 
a particular coordinate direction only leads to a stiffness reduction in the load direction, this approach has 
the potential to more accurately reflect the damage behavior that actually takes place, particularly for 
composites with more complex fiber architectures. 

A wide variety of failure models, which mark the end of the stress-strain curve, have been developed 
for composites. In models such as the Tsai-Wu failure model (Ref. 6), a quadratic function of the 
macroscopic stresses is defined in which the coefficients of the failure function are related to the tensile, 
compressive and shear failure stresses in the various coordinate directions. This model, while 
mathematically simple and easy to implement numerically, assumes that the composite failure surface has 
an ellipsoidal (in 2D) or ovoid (in 3D) shape. In reality, composite failure surfaces often are not in the 
form of simple shapes. More complex models, such as the Hashin model (Ref. 7), also utilize quadratic 
combinations of the macroscopic failure stresses, but utilize only selective terms in the quadratic function 
in order to link the macroscopic stresses to local failure modes such as fiber or matrix failure. However, 
an overall quadratic form to the failure functions (albeit in a piecewise fashion) are still assumed. This 
approach was extended in models such as those developed by Puck et al. (Ref. 8), Pinho et al. (Ref. 9) and 
Maimi et al. (Ref. 10), in which complex equations were developed to predict local failure mechanisms in 
terms of macroscopic level stresses. In this manner, the failure response and complex failure surfaces 
present in actual composites could be more accurately represented. However, in these advanced models, 
very complex tests are often required to characterize the model parameters and the applicability of the 
models may be limited to specific composite architectures with specific failure mechanisms. In a 
combination of approaches, researchers such as Mayes and Hansen (Ref. 11) and Feng (Ref. 12) utilize an 
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approach where stress (or strain) invariants based on macroscopic stresses are used to define the initiation 
of failure. However, different forms of the invariants are used to determine whether fiber-dominated or 
matrix-dominated failure occurs. Given the variety of failure models present in the literature, activities 
such as the World Wide Failure Exercise and its multiple iterations (e.g., Refs. 13 to 15) have attempted 
to conduct a rigorous review of which failure model or models provides the optimum prediction of 
composite failure. In all of these studies, none of the examined models showed a complete ability to 
predict composite failure or displayed a significant advantage compared to the other models examined. 

The difficulty in simulating composite failure can be related to the fact that in reality failure is a 
highly localized phenomenon dependent on various combination of fiber, matrix and interface failures. 
Due to these complex and interacting local failure mechanisms, and the fact that these mechanisms can 
vary based on the constituent materials and fiber architecture, in reality the actual composite failure 
surface often does not conform to a shape that can be easily simulated using a simple mathematical 
function. Conversely, attempts to utilize discrete functions to analyze the complex local mechanisms can 
result in models with a large number of parameters that require a highly complex test program to obtain. 
In the methodology described in this paper, an approach is used in which the actual experimental three- 
dimensional failure envelope for a composite is entered in a tabulated fashion. Specifically, a stress or 
strain based invariant leading to the initiation of failure is defined as a function of the location of the 
current stress state in stress space. In this manner, an arbitrary failure surface can be easily defined based 
on actual experimental data in combination with numerical data obtained using any desired existing 
failure model. The current approach thus serves as a general framework which is not limited to any 
arbitrarily imposed failure surface based on an arbitrarily defined mathematical function. 

In the following sections of this paper, a brief summary of the plasticity-based deformation model is 
presented. The key aspects of the damage model and its characterization are then described. Details of the 
failure model are then discussed, including the overall methodology along with its application for a 
representative composite material. 


Deformation Model Overview 


A complete description of the deformation model is given in Goldberg et al. (Refs. 4 and 5). A 
summary of the key features of the model is presented here. In the deformation model, a general quadratic 
three-dimensional orthotropic yield function based on the Tsai-Wu failure model is specified as follows, 
where 1, 2, and 3 refer to the principal material directions: 


O11 
O22 
fo)=-1+( Fh 0 0 of ® 
912 
23 
Ss (1) 
Fr Fo F3 9 0 90 You 
Fo Fo Fy3 0 0 0 | oy 
+(5}; 522 933 9}. 923 531) eee, tle 
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O° O° O° O Bx. -O Poss 
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In the yield function, 0, represents the stresses and F; and F; are coefficients that evolve and are 
functions of the current values of the stresses in the various coordinate directions. By allowing the 
coefficients to vary, the yield surface evolution and hardening in each of the material directions can be 
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precisely defined. The values of the normal and shear coefficients can be determined by simplifying the 
yield function for the case of unidirectional tensile and compressive loading in each of the coordinate 
directions along with shear tests in each of the shear directions. In the above equation, the stresses are the 
current value of the yield stresses in the normal and shear directions. To determine the values of the off- 
axis coefficients (which are required to capture the stress interaction effects), the results from 45° off-axis 
tests (an arbitrarily chosen angle) in the various coordinate directions can be used. The values of the off- 
diagonal terms in the yield function can also be modified as required in order to ensure that the yield 
surface is convex (Ref. 5). 

A nonassociative flow rule is used to compute the evolution of the components of plastic strain. The 
plastic potential for the flow rule is shown in Equation (2): 


2 2 2 
7 F110} + H 037 + 33033 + 2H 7011079 + 2H 73097033 + 2H3 10330) Q) 
p 2 2 
+ H440{7 + 155033 + 16603} 


where oj; are the current values of the stresses and H; are independent coefficients, which are assumed to 
remain constant. The values of the coefficients are computed based on average plastic Poisson’s ratios 
(Refs. 4 and 5). The plastic potential function in Equation (2) is used in a flow law to compute the 
components of the plastic strain rate, where the usual normality hypothesis from classical plasticity 


(Ref. 16) is assumed to apply and the variable 4 is a scalar plastic multiplier: 


. Oh 
BP =A 3 
eee (3) 


Given the flow law, the principal of the equivalence of plastic work (Ref. 16) can be used to 
determine that the plastic potential function, /, can be defined as the effective stress and the plastic 
multiplier can be defined as the effective plastic strain rate. 

To compute the current value of the yield stresses needed for the yield function, tabulated stress-strain 
curves are used to track the yield stress evolution. The user is required to input twelve stress versus plastic 
strain curves. Specifically, the required curves include uniaxial tension curves in each of the normal 
directions (1,2,3), uniaxial compression curves in each of the normal directions (1,2,3), shear stress-strain 
curves in each of the shear directions (1-2, 2-3, and 3-1), and 45° off-axis tension or compression curves 
in each of the 1-2, 2-3, and 3-1 planes. The 45° curves are required in order to properly capture the stress 
interaction effects. By utilizing tabulated stress-strain curves to track the evolution of the deformation 
response, the experimental stress-strain response of the material can be captured exactly without any 
curve fit approximations. The required stress-strain data can be obtained either from actual experimental 
test results or by appropriate numerical experiments utilizing stand-alone codes. The ability to account for 
rate and temperature effects has also been incorporated into the deformation model. To track the evolution 
of the deformation response along each of the stress-strain curves, the effective plastic strain is chosen to 
be the tracking parameter. Using a numerical procedure based on the radial return method (Ref. 16) in 
combination with an iterative approach, the effective plastic strain is computed for each time/load step. 
The stresses for each of the tabulated input curves corresponding to the current value of the effective 
plastic strain are then used to compute the yield function coefficients. 


Damage Model Overview 


The deformation portion of the material model provides the majority of the capability of the model to 
simulate the nonlinear stress-strain response of the composite. However, in order to capture the changes 
in the unloading modulus observed as the material is unloaded from various stress levels and the local 
softening of the stress-strain response that is often observed in composites (Ref. 17), a complementary 
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damage law is required. Strain equivalence is assumed in the damage law formulation, therefore the total, 
elastic and plastic strains in the actual and effective stress spaces are the same for every time step 
(Ref. 18). The utilization of strain equivalence permits the plasticity and damage calculations to be 
uncoupled, as all of the plasticity computations can take place in the effective stress space. 

The first step in the development of the damage model is to relate the actual stresses to a set of 
effective stresses by use of a damage tensor M: 


o =Mo off (4) 


The effective stress rate tensor can be related to the total and plastic strain rate tensors by use of the 
standard elasto-plastic constitutive equation 


6, =Clé-2,) (5) 


where C is the standard elastic stiffness matrix and the actual total and plastic strain rate tensors are used 
due to the strain equivalence assumption. 

Given the usual assumption that the actual stress tensor and the effective stress tensor are symmetric, 
the actual stresses can be related to the effective stresses in the following manner, where the damage 
tensor M is assumed to have a maximum of 36 independent components: 


O11 on 

O79 on 
eff 

. =[m] 2 (6) 
12 

073 Gat 
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In order to maintain a one to one relationship between the effective stresses and the actual stresses (i.e., to 
ensure that a uniaxial load in the actual stress space does not result in a multiaxial load in the effective 
stress space), the damage tensor is assumed to be diagonal, leading to the following form: 


[Mm, 0 0 0 0 0 
0 My 0 0 0 0 
M]- 0 0 M3 0 0 0 ) 
0 0 0 My 0 0 
0 0 0 0 Ms, 0 
0 0 0 0 0 Me 


An implication of a diagonal damage tensor is that loading the composite in a particular coordinate 
direction only leads to a stiffness reduction in the direction of the load due to the formation of matrix 
cracks perpendicular to the direction of the load. However, as discussed in detail in Goldberg, et al. 
(Ref. 5), in actual composites, particularly those with complex fiber architectures, a load in one 
coordinate direction can lead to stiffness reductions in multiple coordinate directions. To maintain a 
diagonal damage tensor while still allowing for the damage interaction in at least a semicoupled sense, 
each term in the diagonal damage matrix should be a function of the plastic strains in each of the normal 
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and shear coordinate directions, for example the M1; term for the plane stress case has the strain 
dependences shown below: 


Mi, = My f,,2%),eP,) (8) 


Note that plastic strains are chosen as the “tracking parameter” due to the fact that, within the context 
of the developed formulation, the material nonlinearity during loading is simulated by use of a plasticity 
based model. The plastic strains therefore track the current state of load and deformation in the material. 
To explain this concept of damage coupling further, assume a load is applied in the 1-direction to an 
undamaged specimen. The undamaged modulus in the 1-direction is £,, and the undamaged modulus in 
the 2-direction is E,. The stress-strain response of the material is assumed to become nonlinear 
(represented in the current model by the accumulation of plastic strain) and damage is assumed to occur. 
The original specimen is unloaded and reloaded elastically in the 1-direction. Due to the damage, the 
reloaded specimen has a reduced modulus in the 1-direction of pa The reduced area and modulus are a 


function of the damage induced by the loading and resulting nonlinear deformation in the 1|-direction 
(reflected as plastic strain) as follows: 


eft =(-alller e,, (9) 


where d} ‘ is the damage in the |-direction due to a load in the 1-direction. Note that damage is only 
assumed to initiate once plastic strains have started to accumulate. Alternatively, if the damaged specimen 
was reloaded elastically in the 2-direction, due to the assumed damage coupling resulting from the load in 
the 1-direction, the reloaded specimen would have a reduced modulus in the 2-direction of owe The 


reduced modulus is also a function of the damage induced by the load and resulting nonlinear deformation 
in the 1-direction as follows: 


ef! (1a (ep, ea, (10) 


where d a is the damage in the 2-direction due to a load in the 1-direction. Similar arguments can be 


made and equations developed for the situation where the original specimen is loaded in the 2-direction. 

For the case of multiaxial loading, the semicoupled formulation needs to account for the fact that as 
the load is applied in a particular coordinate direction, the loads are acting on damaged areas due to the 
loads in the other coordinate directions, and the load in a particular direction is just adding to the damaged 
area. For example, if one loaded the material in the 2-direction first, the reduced modulus in the 


1-direction would be equal to £ im . If one would then subsequently load the material in the 1-direction, 
the baseline modulus in the 1-direction would not be equal to the original modulus £,,, but instead the 


reduced modulus £4,” . Therefore, the loading in the 1-direction would result in the following further 
reduction in the modulus and effective area (Ajj) in the 1|-direction: 


ef! (ali ep ed?? = (—altler JM abbler ea a) 


aft =(-al} ler, are? =(- aller, J — ables lan 


NASA/TM—2017-219492 6 


These results suggest that the relation between the actual stress and the effective stress should be based on 
a multiplicative combination of the damage terms as opposed to an additive combination of the damage 
terms. For example, for the case of plane stress, the relation between the actual and effective stresses 
could be expressed as follows: 


ll ll 11) eff 

ee=t dil diy dit ost 
22 22 22) _eff 

63 =( aN 2 \t a?) os (12) 
12 12 12) _ eff 

5 =( ai2\i aii di2) 6! 


where for each of the damage terms the subscript indicates the direction of the load which initiates the 
particular increment of damage and the superscript indicates the direction in which the damage takes 
place. For the full three-dimensional case the actual stresses would be functions of damage parameters in 
all six coordinate directions. 

To properly characterize the damage model, an extensive set of test data is required. Due to the 


tabulated nature of the input, each of the damage parameters (d;, , d;, , etc.) has to be determined as a 


function of the plastic strain in a particular coordinate direction (such as ¢/, ). For example, to determine 


the damage terms for the case of loading in the 1-direction, a composite specimen has to be loaded to a 
certain plastic strain level in the 1-direction. The material is then unloaded to a state of zero stress, and 
then reloaded elastically in each of the coordinate directions to determine the reduced modulus of the 
material in each of the coordinate directions. Expressions such as those in Equations (9) and (10) can then 
be used to determine the required damage parameters for the particular value of plastic strain. The process 
needs to be repeated for multiple values of plastic strain in the 1-direction in order to establish the full 
characterization of the variation of the damage parameters as a function of the plastic strain in the 
1-direction. 


Failure Model 


As discussed earlier in this paper, the majority of the available failure models utilize mathematical 
functions to describe the failure surface, which impose a specific shape on the failure surface. An example 
of this concept can be seen in Figure 1. In this figure, a two-dimensional failure surface in the 011-022 
plane for the case of zero shear stresses generated using the Tsai-Wu failure model for a representative 
AS4/3501-6 polymer matrix composite is shown. The properties used to generate the failure surface were 
taken from Daniel and Ishai (Ref. 6) and are listed in Table 1. The two-dimensional version of the Tsai- 
Wu failure model that was used to generate the failure surface is as follows: 


2 2 2 
AiO11 + foS22 + AirS11 + £22922 + S66F12 + 2f1251 1922 =1 (13) 


where oj are the stresses and the coefficients f; and fj; are defined as follows: 


1 1 1 1 

fay Xo I= Yy, fo6 = as 
1 1 1 1 

to “% Yo fo = V¥o fi2 = 5 VSis22 


where X7 is the longitudinal tensile failure stress, Xc is the longitudinal compressive failure stress, Yr is 
the transverse tensile failure stress, Yc is the transverse compressive failure stress and _S is the in-plane 
shear failure stress. Note that the actual failure surface is continuous but slight discontinuities are shown 
in the presented graphs due to numerical issues in the generation of the graph. As can be seen in the 
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figure, the failure surface is elliptical due to the quadratic nature of the equation defining the failure 
surface. In reality, however, as shown for example in some of the cases discussed by Daniel and Ishai 
(Ref. 6), the failure surfaces of actual composites often do not exhibit this simple shape. A schematic of 
this idea is shown in Figure 2, where a representative example of a failure surface of a composite in the 
622-612 plane is shown. This example failure surface cannot be easily defined by a mathematical function 
of the stresses. Alternatively, as shown in the figure, one potential method of defining the points in the 
failure surface is to use a cylindrical type of coordinate system. In this approach, a variable 6 defines the 
relative location of the point on the failure surface in stress space, while a second variable r defines the 
“magnitude” of the failure surface point in the stress space location. Since the relationship between “7” 
and “@’ also cannot be easily defined by a mathematical function for a realistic composite failure surface, 
a tabulated approach, where a series of “7” and “@” pairs are explicitly defined for a given failure surface, 
can provide a more accurate representation of the failure surface. Buyak (Ref. 19) has demonstrated that 
tabulated approaches can be successful in defining complex failure surfaces for advanced materials based 
on experimental data. The tabulated approach allows for the use of experimentally defined failure surface 
data, a failure surface defined using any existing failure model, or a combination of experimental and 
numerically obtained “virtual” data. The combined approach can allow for the case where actual failure 
data is only available for a portion of the total stress space, with “virtual” data being required to fill in the 
gaps in the failure surface. 


Independent and Dependent Variables 


As described above, for the tabulated failure surface definition proposed in this study appropriate 
independent and dependent variables need to be defined. The independent variables need to define the 
location of a point on the failure surface in stress space, and the dependent variable needs to define the 
magnitude of the failure surface point along the lines defined by the independent variables. For the 
current approach, the in-plane and out-of-plane responses will be considered separately. First, the 
definition of the in-plane failure surface will be discussed. For the in-plane failure surface definition, two 
independent variables are defined. The first independent variable will be the ratio of the shear stress to the 
shear failure stress. For selected values of this shear ratio, the location of each defined point on the failure 
surface in stress space is specified by defining the angle of the point in the 01;-022 plane, which becomes 
the second independent variable. This concept is illustrated in Figure | using the Tsai-Wu failure surface 
(for the case of zero shear stress) discussed earlier. Using simple geometric principals, the angle @ for 
each defined point on the failure surface can be defined in terms of the stresses 61; and 022 as shown 
below. To ensure a set of unique, monotonically increasing angles from —180° to 180°, if the stress 622 is 
negative the computed angle is multiplied by —1 as shown below: 


0=cos! tt 5 is 
yO +922 (15) 
Bact =O if on <0 


For the dependent variable, which is used to define the magnitude of the failure surface point given a 
particular location in stress space, a stress invariant first identified by Fleischer (Ref. 20) is used, defined 


as follows for the plane stress case: 
r= 07, +05) +202) (16) 


This invariant can be considered to be like a “radius” from the origin to the failure surface. The factor of 
2 in front of the shear stress term reflects the symmetry of the stress tensor. Stresses or strains can be used 
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in defining the dependent variable, making the model more general. Using an invariant type of term also 
allows for the stress interactions to be more appropriately accounted for in the failure definition and helps 
to ensure that the failure definition will be accurate for a variety of loading conditions. The general 
concept of the “radius” of the failure surface is shown in Figure 1. 

One difficulty in using a tabulated approach of this type is that for many polymer matrix composites 
the failure surface is highly anisotropic due to the high failure strength in the fiber direction and the much 
lower failure stresses in the transverse direction. This anisotropy can be observed graphically in Figure 3, 
in which a version of the failure surface shown in Figure 1 is shown. In this plot, the range of the y axis is 
adjusted to be equal to the range of the x axis. As can be seen in this figure, the actual failure surface is 
highly elongated. The consequences of the significant elongation of the failure surface can be observed in 
Figure 4, where the radius versus the angle is plotted for the failure surface defined above. As can be seen 
in the figure, the majority of the points of the failure surface are clustered at angles near —180°, 0°, and 
180°. This clustering would create round-off and interpolation difficulties in the numerical 
implementation of the failure model, making the current definition undesirable. 

To facilitate a more even distribution of the tabulated failure surface points along the entire spectrum 
of angles, the stresses used in the definition of the independent variable 6 can be scaled. For the current 
method, the stresses in the 1 1-direction are arbitrarily scaled by the longitudinal tensile failure stress, and 
the 22-direction stresses are scaled by the transverse tensile failure stress. Revised angles are then 
computed using these scaled stresses as shown below. The value of the dependent variable r, however, is 
still computed using the actual unscaled stresses: 


_ SO 
©11-scaled = xX 
T 
O22 
©22-scaled = Y. 
° (17) 
G= cos. ©} |~scaled 


2 2, 
\o; 1-scaled + 922-~scaled 
Bact = Vif 555 <0 


The modified version of the failure surface defined in Figure 1 is shown in Figure 5. The overall shape of 
the failure surface is identical, only the relative values of the axes have changed. A revised plot of the 
radius versus the scaled angle is shown in Figure 6. As can be seen in this figure as compared to Figure 4, 
the points defining the failure surface are much more evenly distributed along the range of angles and the 
sharp points observed in the graph in Figure 4 are smoothed out. This relationship, particularly when 
transferred to a tabulated format, will allow for much smoother and more accurate interpolations when 
used in a numerical implementation. 

As mentioned earlier, for the proposed approach the procedure described above needs to be repeated 
for several ratios of the shear stress to the shear failure stress in order to fully define the in-plane failure 
surface. An example of the scaled failure surface in the 611-022 plane where the shear stress is equal to 
50 percent of the shear failure stress is shown in Figure 7. In this plot, the failure surface for the case with 
zero shear stresses is shown in green for comparison. In this plot, one can observe that increasing the 
shear stresses results in a failure surface with the same shape as the surface generated assuming zero shear 
stresses, but the location of the failure surface is shifted somewhat towards the negative stress quadrant of 
the graph. A plot of the radius versus the scaled angle for the 50 percent shear case is shown in Figure 8, 
along with the plot for the zero shear stress case for comparison (shown in green). Comparing the two 
plots, the relative shapes of the two curves are similar, but the magnitudes of the radius values for the 
50 percent shear stress case are somewhat lower. 
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Modified Center of Failure Surface 


As can be observed in Figure 7, the failure surfaces are not centered at the origin, and as the shear 
stress is increased the failure surface becomes even more skewed away from the origin. At high values of 
shear stresses, the failure surface may not even include the origin. Even when the origin is included, the 
fact that the failure surface is not centered on the origin may cause the angle calculations to be skewed. If 
the failure surface does not include the origin, the angle calculations would not even be valid as they 
would not be unique. 

To mitigate this issue, in the current approach the origin of the scaled failure surfaces for the selected 
shear stress ratios is redefined such that it lies in the center of the failure surface. The modified center of 
the failure surface is defined based on the maximum and minimum scaled stress values in the 11- and 
22-directions as shown below. Modified stresses in the 11- and 22-directions are then defined in terms of 
the modified center as shown below: 


1 
©} 1~center = ~(o, I-max + 9] ey 
2 
1 
©22-center = 5 (©22-max + 599- min) (18) 


0711—mod = ©11-scaled ~ ©11-center 


022-—mod = © 22-scaled ~ 9 22-center 


The angle calculations shown in Equations (15) and (17) are carried out using the revised stresses. The 
radius calculations are carried out using the original, unscaled and unmodified stresses. The modified, 
scaled, re-centered failure surfaces are shown in Figure 9, where the modified surface for the case of zero 
shear stress is shown in green and the modified failure surface for the case where the shear stresses are 
equal to 50 percent of the shear failure stress is shown in black. In both cases, the overall shape of the 
failure surface does not change from the original case but the failure surfaces are now centered about the 
modified origin. 

The radius versus scaled angle plots for the zero shear case (green curve) and for the 50 percent shear 
ratio case (black curve) determined using the modified failure surface centers are shown in Figure 10. In 
both cases, the tabulated points defining the failure surface are more evenly distributed compared to the 
cases with the original origin location. This improvement in the failure point distribution will improve the 
implementation of the failure model. 


Algorithm for Implementation of In-Plane Failure Model 


To apply the in-plane version of the failure model described above, first the failure surface for a given 
composite needs to be converted into a tabulated form in order to establish a benchmark for determination 
if failure has occurred for a given load condition. The baseline failure surface for a given material can be 
established based on experimental data, any available numerical failure model, or a combination of the 
two. In the combination approach, experimental data can be used where available, and analytically 
generated data can be used to fill in the required gaps. Given the baseline failure surface data, 

Equations (15) to (18) are then used to generate a series of “radius” versus “angle” plots for various ratios 
of shear stress to shear failure stress, similar to Figures 13 and 14. For use in a numerical algorithm, the 
data needs to be transformed into a set of tables which serve as input to the computer model. The stress 
values defining the “center” of the failure surface for the selected shear stress ratios also need to be 
defined for use in the calculations. 
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The specific algorithm to implement the failure model is as follows: 


Given: applied load condition 611, 622, 612 
Compute the shear ratio R = 612/S 
3. Scale the longitudinal and transverse stresses: 


NO Re 


O11 
011-scaled = Xx. 
T 
(19) 
O22 
©22-scaled = Y. 
T 


4. Look up the appropriate modified failure surface center stresses (011-center ANd 022-center) 
(Eq. (18)) corresponding to the current shear ratio R. Note that interpolation of these stresses may be 
required if the current shear ratio R is not explicitly defined in the input tables. 

5. Compute the modified longitudinal and transverse stresses given the stresses defining the modified 
failure surface center: 


O11—-mod = O11-scaled ~ 91 1—center (20) 


092-mod = 922-scaled ~ © 22-center 


6. Compute the angle 0 associated with the modified longitudinal and transverse stresses: 


Haas! O11—mod 
2 2 21 
O11—mod + 922-mod ey 
Bact =-0if ©27-mod © 0 
7. Compute the radius associated with the unscaled, unmodified applied stresses: 
r= yo}, +03, +264) (22) 


8. Compare the radius computed at step 7 to the tabulated radius at failure for the given shear ratio R and 
angle @. Note that interpolation of the radius at failure may be required if the specific values of R and 
6 are not specified in the input tables. 

9. Ifthe radius computed using the applied stresses is equal to or greater than the radius at failure, the 
material is assumed to have failed. In the numerical implementation, the element will either be eroded 
once failure is deemed to have occurred or have all of its stiffness properties reduced to negligible 
values. A progressive reduction in stresses will not occur. 


Out-of-Plane Failure Surface Definition and Modification of Failure Model 


The methodology and procedures described above are valid for the case of in-plane loading and 
stresses. For thin composites, particularly under in-plane loading, this case is most likely sufficient since 
the out-of-plane normal and transverse shear stresses can be assumed to be negligible. However, for 
thicker composites, particularly under impact loading conditions, the out-of-plane stresses can often be 
significant and must be included. One way to account for failure due to out-of-plane stresses is to utilize a 
delamination model, as is discussed in many sources such as Barbero (Ref. 17). In this case, the in-plane 
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failure model could be used in combination with a delamination model to fully define the composite 
failure. However, an alternative approach would be to define a tabulated failure surface specified by the 
out-of-plane stresses that is independent from the failure surface defined by the in-plane stresses. As will 
be described below, the out-of-plane failure surface can be defined in a manner similar to what was used 
to define the in-plane failure surface. The radius-like dependent variable can then be used to couple the 
failure responses due to the in-plane and out-of-plane stresses. 

To generate the required out-of-plane failure surfaces, a series of two-dimensional failure surfaces in 
the 613-023 plane are defined for various ratios of the out-of-plane normal stress 033 to the out-of-plane 
normal failure stress. Equations similar to Equations (15), (17), and (18) can then be used to define the 
values of the angle 6 (where o1; in the equations is replaced with o13 and 622 in the equations is replaced 
with 623) for the out-of-plane failure surfaces. For the dependent variable, a stress invariant radius term 
similar to Equation (16) is defined: 


r= 02, +207, +203 23 
33 13 23 


where the transverse shear stresses are multiplied by two due to the symmetry of the stress tensor. Plots of 
the radius versus scaled angle (similar to those shown in Fig. 10) can then be generated for several out-of- 
plane normal stress ratios based on the out-of-plane failure stresses and converted into a tabulated set of 
input for the failure model. For a given out-of-plane loading condition, the out-of-plane normal stress to 
out-of-plane normal failure stress ratio, angle and radius can be computed using a process similar to the 
process used for the in-plane stresses. The radius computed using the current stresses can then be 
compared to the corresponding radius at failure to determine if failure has occurred. 

For the case where both in-plane and out-of-plane stresses are present in a material under a specified 
load state, a method needs to be specified to allow for the stress interactions that may take place under the 
combination of in-plane and out-of-plane loading. Accounting for the interaction is important because the 
combination of stresses may cause failure in the material even though the individual in-plane and out-of- 
plane criteria may not indicate that failure has occurred. To account for the stress interaction, a ratio dip 
and a ratio doo» are defined. The ratio d;, represents the ratio of the radius computed using Equation (16) 
for the applied in-plane stresses to the corresponding radius computed using the in-plane stresses at 
failure for a given angle and shear stress ratio. The ratio d,op represents the ratio of the radius computed 
using Equation (23) for the applied out-of-plane stresses to the radius at failure for a given out-of-plane 
angle and out-of-plane normal stress ratio. A weighted sum of the ratios can then be computed using the 
following expression, where n is an arbitrary number defined by the user: 


d="d" +d", (24) 


If the weighted ratio d is greater than or equal to one, failure is assumed to have occurred. If the weighted 
ratio d is less than one, failure is assumed to have not occurred. 


Conclusions 


A generalized composite model suitable for use in polymer composite impact simulations has been 
developed. The model utilizes a plasticity based deformation model obtained by generalizing the Tsai-Wu 
failure criteria. A strain equivalent damage model has also been developed in which loading the material 
in a particular loading direction can lead to damage in multiple coordinate directions. A general, tabulated 
failure model has also been formulated. Within the failure model, a methodology has been developed to 
convert two-dimensional failure surfaces to a tabulated format based on the location of the points of the 
failure surface in stress space and a dependent stress invariant variable which defines the magnitude of the 
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failure stress points at a given location. An algorithm to implement the failure model has been developed, 
and a procedure for accounting for a combination of in-plane and out-of-plane stresses has been specified. 

Future efforts will include developing the detailed numerical algorithms required to implement the 
damage and failure models into the material model being developed for inclusion within the LS-DYNA 
computer code. The failure model will also be refined to include appropriate techniques for element 
removal once the failure criteria has been satisfied. Methods to account for delamination failure will also 
be formulated. Extensive sets of verification and validation studies will also be undertaken in order to 
fully exercise the developed model. 
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TABLE 1.—MATERIAL PROPERTIES FOR AS4/3501-6 
UNIDIRECTIONAL COMPOSITE 


Property Value 
Longitudinal Modulus £1 (GPa) 147 
Transverse Modulus E22 (GPa) 10.3 
Longitudinal Poisson’s Ratio viz 0.27 
Transverse Poisson’s Ratio v23 0.54 
In-Plane Shear Modulus Gi2 (GPa) 7.0 
Longitudinal Tensile Strength X7 (MPa) 2280 
Longitudinal Compressive Strength Xc (MPa) 1725 
Transverse Tensile Strength Yr (MPa) 57 
Transverse Compressive Strength Yc (MPa) 228 
In-Plane Shear Strength S (MPa) 76 
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Figure 1.—Tsai-Wu failure surface for AS4/3501-6 composite. Variables @ and r used to specify tabulated 
failure surface are also defined. 
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Figure 2.—Schematic of failure surface for representative composite. 
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Figure 3.—Tsai-Wu surface for AS4/3501-6 composite with equal x and y axis ranges. 
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Figure 4.—Plot of radius vs. angle for original AS4/3501-6 failure surface. 
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Figure 5.—Scaled Tsai-Wu failure surface for AS4/3501-6 composite. 
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Figure 6.—Plot of radius versus angle for scaled AS4/3501-6 failure surface. 
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Figure 7.—Scaled Tsai-Wu failure surface with shear stress equal to 50 percent of 
shear failure stress for AS4/3501-6 composite (black curve). Failure surface with 
zero shear stress shown in green. 
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Figure 8.—Plot of radius versus angle for scaled AS4/3501-6 failure surface for case 
with shear stress equal to 50 percent of shear failure stress (black curve) and for 


case of zero shear stress (green curve). 
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Figure 9.—Scaled Tsai-Wu failure surface with modified center for AS4/3501-6 composite. 
Surface for case with zero shear stress shown in green, surface for case with shear stress 
equal to 50 percent of shear failure stress shown in black. 
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Figure 10.—Plot of radius versus angle for scaled failure surface with modified center 
for AS4/3501-6 composite for case of zero shear stresses (green curve) and case 
where shear stresses equal 50 percent of shear failure stress (black curve). 
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